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Abstract: A diagonalization scheme for the Rabi Hamiltonian, which describes a qubit 

(as; 

interacting with a single-mode radiation field via a dipole interaction, is proposed. It is 



shown that the Rabi Hamiltonian can be solved almost exactly using a progressive scheme 
that involves a finite set of one variable polynomial equations. The scheme is especially 



1 efficient for the lower part of the spectrum. Some low-lying energy levels of the model with 
Oh! 

several sets of parameters are calculated and compared to those provided by the recently 
C ■ 

d 1 proposed generalized rotating-wave approximation and a full matrix diagonalization. 

3 . 
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It is well known that a qubit (two-level atom) interacting with a radiation field is a 



O 

OO 
(N 

. fundamental problem in modern physics, encountered, for example, in condensed matter, 
biophysics, and quantum optics. I 1 ' 2 ' The most common approximation assumes that the ra- 
diation field is quasi-monochromatic, which leads to the so-called Rabi Hamiltonian. t 3 '4 Due 
to its simplicity and importance, the Rabi Hamiltonian has been studied extensively using 
various methods. More complicated extensions to many two-level atoms interacting with 
a quasi-mono-chromatic field are described by the Dicke Model, t 5 ' The Rabi Hamiltonian 
is also equivalent to that of the spin-boson model, ' 2 ' 6 ' which describes a qubit interacting 
with an infinite collection of harmonic oscillators that model the environment acting as 
a dissipative bosonic bath. The latter model has also attracted considerable attention in 
various quantum many-body systems due to the rich physics of quantum criticality and 
decoherence.t 7-9 ' Though it has long been conjectured that the Rabi Hamiltonian may 
be exactly diagonalizablej 1 ^ with many analytical and numerical studies carried out for 
decades, I 11-21 ] exact solutions have not been presented except for some special cases J 22-31 ! 
Due to the absence of an exact treatment, approximations have been developed. The rotating 
wave approximation (RWA),! 32 ] which neglects the counter- rotating term in the Hamiltonian, 
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a case with the counter-rotating term, and another case with degenerate level of the atom 
are three well-known solvable cases. [ 30 > 31 1 

The Rabi Hamiltonian can be written as' 3,4 ' 32 ! 



Hi 



cu a^a + ujS + g(a i + a)(S- + S+), 



(1) 



where uq and u are frequencies of a single-field mode and the level splitting of the two- 
level atom, respectively, g is the coupling parameter between the atom and the field, a 
and are annihilation and creation operators of the field mode, respectively, satisfying the 



usual bosonic commutation relation: [a, a'] = 1, [a, a] = [ 



a\ a' 



0, So is the operator of 



atomic inversion, and S + (S-) is the atomic raising (lowering) operator, which satisfy the 
SU(2) Lie algebraic relations [So, S±] = ±S±, [S + , S-] = 2S - The pseudo-spin operator 
S 2 = S + S- + So(So — 1) and the parity P = exp (i-ir(a^a + So + S)\ where S — 1/2 is 
the pseudo spin of the atom, commute with the Hamiltonian (1), namely [S 2 ,H-r] = and 
[P,Hr] = 0. 

Let |0;/i) (fi = + or — ) be simultaneously the boson vacuum and the eigenstate of So 
satisfying a\0; fi) = 0, So|0;±) = ±^|0;±). The Hamiltonian (1) under the basis spanned 
by {\k; fj) = -^ya^ fe |0; fi); k — 0, 1, 2, • • ■} has the following tridiagonal form: 



H 



R 



(2) 





£T(-) 

of which the order of the Fock states is |0; +), |1; — ), \2; +), • • •, • • •, |0; — ), |1; +), \2; —) 

■ ■ •, where the two tridiagonal infinite sub-matrices and are given by 



f±\u gVl ••• 

gVl w t|w g^/2 ••• 

gV2 2uo±\u gV% 0--- 

V ; ••• ••• ••• ••• J 



(3) 



Under the afore arranged order of the basis, the effective Hamiltonian (3) may be written 



as 



= (ujQcfia + g(a + a^)) ± —oj coa{rca)a). 



(4) 



The spin component fi = + or — in the Fock states \k; fi) will be omitted in the following 
for simplicity. The first term in (4) is easily diagonalizable under the shifted boson states 
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| n) = e -\{ 9 /^? y^( Q t + g / UQ )" e -% a * | ) (5) 

with the corresponding eigenvalues E n = uqu — g 2 /ojq for n — 0, 1, 2, ■ • -. The second term 

is diagonal in the original boson subspace spanned by the Fock states Therefore, the 
Hamiltonian (4) can also be expressed as 



OO ,,00 

UJ 



# (±) = £ E n \n){n\ ± — ^fu^2\2k + I) {2k + 1|. (6) 

n=0 ^ k=0 

Since the structure of the Hamiltonians and fl^ - ) is the same, in the following, we 
take to elucidate our diagonalization procedure, and assume that the parameters uo, 
u, and g are all non-zero. We adopt the analytical step-by-step diagonalization procedure 
proposed in [33]. In the k-th. step, we take a diagonalized part of the Hamiltonian H^ -1 ' in 
the shifted boson basis (5) and the k-th projection in the second term of (6) with H^ k ' = 
jjik- 1 ) — uj\2k + l)(2k + 1| to do the diagonalization. In the initial k = step, 

H^ = J2E n \n)(n\ + ^-u\l)(ll (7) 

n=0 Z 

which can be diagonalized exactly with the eigenstates 



a 1 



-WE — — (8) 

\ N T0 n=0 E n + uj/2 — E T0 
where is the normalization constant, a 2 l k+1 is the overlap of the original Fock state 
\2k + 1) with the shifted state \n) given by 



v 2 n k+l = (n\2k + 1) = V / (2A; + l)!n!e- ( 4 )2 (A)«-2^i £ 



2q 



W Y( n ~ 2A; - 1 + 9) ! (2A;+l-g)!g!' 

(9) 

One can directly check with the eigenequation H^°'\^^) = E^\ip!^) that (8) is indeed 
the eigenstate of with the corresponding eigenvalue E^ if the following equation is 
satisfied: 



fa 1 ^ 2 



i-^e - v * „ (0) =°» ( io ) 



n=0 £?„ + w/2 - £• 



where To labels the ro-th root of Eq. (10). However, as shown in Fig. 1, the overlaps ce%° +1 
with fixed k are only non- negligible for a finite number of values n with n = mi, mi + 
m 2 . As a result, only a finite number of states {\n)} correlate among each other in 
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the k-th step, which, in turn, effectively truncates the infinite sum in (10) into a finite sum 
of those terms with non-negligible a^ k+1 for n = mi, mi + !,■■■, m 2 . The number of values 
n defined by dk = m 2 — mi + 1 is called correlation length. While other states with n > m2 
or < mi remain essentially unchanged after the diagonalization. The correlation length dk 
increases with increasing k and the coupling parameter g/ojQ. Typically, dk ~ 10-20 when 
g/ujo < 1, and dk ~ 20-60 when 1 < g/uo < 2. Furthermore, only those states {\n)} within 
the correlated region will be correlated among each other. In short, the matrix form of 
Hamiltonian in the shifted boson basis is almost block diagonal. As a consequence, Eq. 
(10) can effectively be written as 



1 



do-l 

-E 

n=0 



a 



1\2 



E n + oj/2 



-C/tq 



0. 



(11) 



where mi = and m 2 — do — 1. After the initial step, given by (6) can be rewritten as 



'do-l 



H M = I E E^)(^\ + E ^W(*l I 

^to=0 n=do+l / fc=l 



]T |2fc+l)(2Jfe + l| 



;i2) 



(a) 



(b) 



0.1 r 
0.08 - 
0.06 
0.04 
0.02 
0. - 



0. — 

10 15 20 10 20 30 40 50 60 
0.15 1 i 



(0 



0.05 



(d) 



20 40 



10 20 30 40 50 60 



Fig. 1. Overlaps functions of n for several values of k and the coupling parameters g/ujQ. 

(a) k = 2, g/uo = 0.5; (b) k = 12, g/u Q = 1.0; (c) k = 15, g/oj = 2.0; (d) k = 3, g/wo = 2.5. It 
can be seen that the corresponding correlation length is (a) d 2 ~ 10, (b) di 2 ~ 30, (c) di§ ~ 50, 
and (d) di ~ 25, respectively. 
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where 



do-l 

E 



a, 



\ N T ° ] fo E n + u/2 - E% ] 



\n). 



(13) 



Next, using a similar procedure, the first term and the k = 1 projection in the second 
term of (12) given by 



do-l 



H {1) = E EiX)^ + E E n \n){n\ -u)\3)(3\ 



TO I T TO 

ro=0 n=do 



(14) 



needs to be diagonalized. If the correlation length d\ > d , the correlated eigenstates may 
be written as 



1 



\ (1) 

\ Nti \tq=1 ^to — Cjt\ n=d 

where it is assumed that (i(;^\3} ~ when r < I, and 



TO=l ^TO 



tin n=do &n + U/Z — hi- 



(1) 
Tl 



\n 



(15) 



(^ } |3) 



i do-l 

E 



1 ^ 

ol or 



(16) 



according to (13), while the uncorrelated eigenstates satisfy |"0t ) = with the corre- 

sponding eigenvalue E^ for r < I, and = \n) with the corresponding eigenvalue E n 

for n > d\. Using (14) and (15) to solve the eigenequation H^\%jj^) = E^\ip^), one can 
prove that (15) is indeed the eigenstate of if the corresponding eigenvalue E^ is the 
7~i-th root of the following equation: 



1 ^ ~^(o) — Tai) 2- 



a 



3\2 



3 =Z Erg' — Erf n=d E n + 00 / 2 — E T1 



(1) 



0. 



(17) 



According to this procedure, once eigenvectors IV^) and the corresponding eigenvalues 
E^ are known from the k-th step, the k + 1-th step results are given by 



-l) E 



U» {^\2k + l) 



(18) 



for Tfc + i = mx,mi + 1, • • • , m 2 and c4+i = m 2 — mi + 1, in which should satisfy 

™ ((^|2A;+1)) 2 n 
l-oo ^ 



(19) 



While — ) w ith the corresponding eigenvalue = for 7fc+i = 

0, 1, • • ■ , m\ — 1, and r^+i = m 2 + l,m 2 + 2, • ■ -. In fact, with increasing k, the region 
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TABLE I: The first few low-lying energy levels (ordered from the lowest) determined by the Rabi 
Hamiltonian (1) obtained from the Progressive Diagonalization Scheme (PDS) are compared to 
those calculated using the GRWA and results from a Full Matrix Diagonalization (FMD), which is 
done by diagonalizing the infinite tridiagonal matrix (3) in a truncated Hilbert subspace. 





= 0.5, UJQ 


= 0J 


g/u 


= 1.0, uj 


= u) 


g/^o = 


2.0, w = 


= 0.75w 


PDS 


GRWA 


FMD 


PDS 


GRWA 


FMD 


PDS 


GRWA 


FMD 


-0.6332 


-0.5532 


-0.6333 


-1.1479 


-1.0677 


-1.1479 


-3.0225 


-3.0002 


-3.0225 


-0.1200 


-0.0602 


-0.1200 


-1.0100 


-0.9936 


-1.0102 


-3.0221 


-2.9999 


-3.0222 


0.6953 


0.7401 


0.6954 


-0.2317 


-0.2110 


-0.2317 


-2.2793 


-2.2546 


-2.2794 


0.8253 


0.8635 


0.8253 


0.1334 


0.2642 


0.1334 


-2.2738 


-2.2474 


-2.2739 


1.5870 


1.5977 


1.5870 


0.9270 


0.9403 


0.9270 


-1.5524 


-1.5142 


-1.5525 


1.9355 


1.9115 


1.9355 


1.1048 


1.0539 


1.1048 


-1.5172 


-1.4966 


-1.5172 


2.5485 


2.5375 


2.5485 


1.8428 


1.8559 


1.8429 


-0.8554 


-0.8385 


-0.8554 


2.9477 


2.9528 


2.9478 


2.1436 


2.1269 


2.1437 


-0.7376 


-0.6792 


-0.7376 



of the correlation among {|n)} moves to larger values of n, and keeps lower part of the 
levels with small n less affected. Therefore, after a number of steps, several of the low-lying 
eigenstates of the Hamiltonian (6) and the corresponding eigenvalues can be obtained. The 
results are almost exact because the overlaps (ip^\2k + 1) are small enough to be neglected 
when Tfc satisfies (ifj^\2k + 1} ~ 0, which is similar to the situations shown in Fig. 1 for 
a n fc+1 = {n\2k+l) because the overlaps (ip^\2k+l) is a linear function of = {n\2k+\). 

Since the Rabi Hamiltonian can not be solved exactly except the case with the RWA [32] 
and other two special cases studied in [30-31]. As shown in [20], the Generalized Rotating- 
Wave Approximation (GRWA) is significantly more accurate than the RWA. To demonstrate 
our Progressive Diagonalization Scheme (PDS), a few low-lying energy levels with three sets 
of model parameters are calculated and compared to those calculated using the GRWA and a 
Full Matrix Diagonalization (FMD), which is done by diagonalizing the infinite tridiagonal 
matrix (3) in a truncated Hilbert subspace. The dimension of the truncated subspace, 
with which the first few low-lying energies keep unchanged with further increasing of the 
dimension, depends on the coupling parameter g/ojQ. For g/u>o = 0.5-2.0, the dimension of 
the truncated subspace is 10 3 — 10 5 , with which the FMD can be preformed on a PC. When 
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g/u>o > 2.0, in order to get accurate results, the FMD needs more computer memory to store 
matrix elements and many hours of CPU time to do the diagonalization. Therefore, we only 
compare the PDS with the GRWA and the FMD for g/uo < 2.0. The results are shown 
in Table 1. The PDS results were obtained with less then 4 steps for g/uo = 0.5, 14 steps 
for g/uo = 1.0, and less than 30 steps for g/uo = 2.0, which are exactly the same as the 
corresponding results obtained from a full matrix diagonalization. In each step, only a single 
variable finite order polynomial equation needs to be solved. The order of the polynomial 
equations used in the PDS depends on the correlation length df.. In the cases shown in 
Table 1, the maximal dk is around 60 when g/uo = 2.0. Such polynomial equations can 
easily be solved by using Mathematica with a few seconds of CPU time. Though the PDS 
for g/uo = 2.0 needs about 30 steps in order to get accurate results, the total CPU time 
needed is much less than that of the FMD. Besides the FMD needs more memory, the CPU 
time needed is 10-20 times more than the total CPU time needed for the PDS. The PDS 
results are certainly better than those obtained from the GRWA, which nonetheless has the 
advantage of being a simpler to implement. 

In conclusion, a diagonalization scheme of the Rabi Hamiltonian in the shifted boson 
basis is proposed based on the step-by step diagonalization method, from which lower part 
of eigenstates and the corresponding eigenvalues can be obtained progressively within a finite 
number of steps. In each step, only a single variable finite order polynomial equation needs to 
be solved due to the fact that only a finite number of the shifted boson states correlate among 
each other. Such type of polynomial equations also appear in the Tamm-Dancoff and random 
phase approximations with separable potentials in nuclear physics. ^ The convergence is 
related to the number of steps that are needed, which in turn depends on the coupling 
parameter g/uo and and position of the energy levels to be calculated. Only a few steps are 
needed for g/uo < 1, while 10-30 steps are needed in order to get accurate results for the 
low- lying part of the spectrum for 1 < g/ujQ < 2. More steps are needed for higher excited 
states. Anyway, the total CPU time needed is much less than that of the truncated matrix 
diagonalization for the low- lying part of the spectrum. As shown in [20], the coupling g/cjQ 
is often small with g/ujQ < 1 in quantum optics, microwave resonator, superconducting LC 
circuits, but g/ooo ~ 1 for nano-mechanical resonator coupled to a charged qubit. Higher 
values of g/uo may also be possible for spin-boson models of other systems. Since the PDS 
scheme yields nearly exact results that seem to be independent of the model parameters Uq, g, 
and uj, and the single variable finite order polynomial equation involved in each step is much 
more easily to be solved than the truncated matrix diagonalization, it should be a useful 
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scheme for studying dynamics in the spin-boson model for strong coupling cases, especially 
in studying quantum critical phenomena and decoherence, where accurate solutions to the 
problem is essential. Related work is in progress. 
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